function[SPMatrix,ttf,freq1,C,Coh,TF,C1,Err] = induction(),

%---------Pl do not change the frequency setting !!!!-------------------
%TLFreq = [1 2 3 4 5];
%TLFreq = [0.4668 1.1288    1.2743    1.4384    1.6238    1.8330    2.0691 2.3357    2.6367    2.9764    3.3598    3.7927    4.2813    4.8329 5.4556 ];
%TLFreq = [0.05 0.1 0.2 0.2825    0.5531    1.3966    1.8979    2.1684    2.7891 3.2825    3.6167    3.8395    4.2294    4.4125    4.6114    5.2162 5.7175];
%TLFreq = [0.1000    0.1565    0.2448    0.3831    0.5995    0.9380    1.4678   2.2967    3.5938    5.6234]; %logspace(-1,0.75,10)
%TLFreq = [0.1000  0.15 0.3 0.7278    1.3556    1.9833    2.6111    3.2389    3.8667  4.4944  5.1222 5.8];
%TLFreq = logspace(-1,0.75,15);
%TLFreq = linspace(2,5.5,15);
%TLFreq = [0.01 0.02 0.03 0.04 0.05 0.08 0.1000  0.15 0.3];
%TLFreq = logspace(-3,-0.515,10);
%TLFreq = [0.0010  0.0036  0.0127    0.0240    0.0454  0.0857    0.1618    0.3055];
%TLFreq =  [0.0010    0.0016    0.0025    0.0039    0.0062    0.0097    0.0153    0.0241    0.0379    0.0597    0.0941    0.1483    0.2336    0.3680   0.5796    0.9131    1.4384    2.2660    3.5697    5.6234];
%---------Pl do not change the frequency setting !!!!-------------------


%-------------2 hours data
% data = load('c:\manoj\projects\induction\Test_01_02h.txt');
% ProcDef.block = 512;
% ProcDef.overlap = 1;
% Parzen = 0.5;
% DeltaFreq = [12/1024];%12 = 24/2; 1024 as a standard fft length
% flag = 0; %flag = 0 for robust band flag = 1 for parzen
% trash = 0.3; %Coherency trashold (between 0 and 1). 0 means mean
% %TLFreq = [0.1000  0.15 0.3 0.7278    1.3556    1.9833    2.6111    3.2389    3.8667  4.4944   5.8];
% TLFreq = [0.0941    0.1483    0.2336    0.3680   0.5796    0.9131    1.4384    2.2660    3.5697    5.6234];
%---------------------------------------------

%--------------4 hours data
data = load('c:\manoj\projects\induction\Test_01_04h.txt');
ProcDef.block = 2048;
ProcDef.overlap = 1;
Parzen = 0.6;
DeltaFreq = [6/1024]; % 6 = 24/4; 1024 as a standard fft length
flag = 0; % flag = 0 for robust band flag = 1 for parzen
trash = 0.8;% Coherency trashold (between 0 and 1). 0 means mean
%TLFreq = [0.0010  0.0036  0.0127    0.0240    0.0454  0.0857    0.1618    0.3055];
TLFreq = [0.0039    0.0062    0.0097    0.0153    0.0241    0.0379    0.0597    0.0941    0.1483    0.2336    0.3680];
%---------------------------------------------------

nchannel = 2;
len = floor(length(data)/ProcDef.block);
ProcDef.SelStacks = ones([1,len]);
[TIndex,len] = GetSelArray(ProcDef);
ProcDef.len=len;
ProcDef.SelStacks = ones([1,len]);
[TLine,TLRad,TLFreq] =  SetTLine_S(TLFreq,DeltaFreq,Parzen,ProcDef.block);
nfrq = length(TLFreq);
ProcDef.Weights = ones([nfrq,len]);
disp('GetRobSpectra -> Now making matrices of cross and auto spectra....');
WI = hanning(ProcDef.block);
freq1 = TLFreq;

SPMatrix = zeros([nfrq len nchannel nchannel]);

for i = 1:len,

for u = 1:nchannel,
   K = CorTre(data(TIndex(i,1):TIndex(i,2),u),0);
	Channel(:,u) = fft(K.*WI);
end; 	
% auto and cross spectra---------------------------
for ii = 1:length(TLFreq),
   par = parzen1(TLRad(ii)); % Bit inefficient (as it is called many times to do the same job)
   SPMatrix(ii,i,:,:) = CompuSPM1(Channel,par,TLine(ii),TLRad(ii),flag);
end;
%------------------------------------------------
%clear Channel;

end;


TF = ind_tf(SPMatrix);
Coh = cohr(SPMatrix);
LL = Coh > trash;
a = size(SPMatrix);
if a(2) > 1,
    for ip = 1:length(TLFreq),
    SPM(ip,:,:) = mean(squeeze(SPMatrix(ip,LL(ip,:),:,:)),1);
    end;
else,
    SPM = squeeze(SPMatrix);
end;

%Variance

for ik = 1:length(TLFreq),
    df = TLRad(ik)*2*sum(LL(ik,:)); %number of independent frequencies
    if df < 5,
        df = 5;
    end;
    fis = fischer(df-4); % calculate fischer distribution F (4,df-4,0.05) Bendat & Piersol, page 110
    Err(ik) =  (4/(df-4)) * fis * (1.0 - mean(Coh(ik,LL(ik,:))) )* SPM(ik,1,1)/SPM(ik,2,2); % Power of output
    % There is a problem here. The mean coherency of the SPCTRA is computed not the
    % coherency of the mean SPECTRA !
end;
ttf = sum(TF.*LL,2)./sum(LL,2);
TF = ind_tf(SPM);
a = 6371.2;    %a = earth surface
L = 1; %Spherical harmonics degree = 1
C = a*(L - (L+1)*ttf)./(L*(L+1)*(1+ttf)); %Weidelt's admittance
C1 = a*(L - (L+1)*TF)./(L*(L+1)*(1+TF)); %Weidelt's admittance
ErrC = a*(L - (L+1)*Err)./(L*(L+1)*(1+Err));
disptf(C,ErrC/10,freq1,'k');


function[SPM] = CompuSPM1(Channel, par,TL,TR,flag),

a = size(Channel);
nchannel = a(2);

for i = 1:nchannel,
   for j = 1:i,
       if flag == 1,
           talkal = sum(Channel(TL-(TR-1):TL+(TR-1),j).*conj(Channel(TL-(TR-1):TL+(TR-1),i)).*par); % summing w.r.t parzen window%
       elseif flag == 0,
           talkal = rbpz(Channel(TL-(TR-1):TL+(TR-1),j).*conj(Channel(TL-(TR-1):TL+(TR-1),i))); % Robust Band averaging
       end;
      if i ~= j,
       SPM(i,j) = talkal; %Note the real and imaginary are
       SPM(j,i) = talkal; %note sep. packed like MT processing !
      elseif i == j,
       SPM(i,j) = talkal;
     end;
  end;
end;

function[PF] = parzen1(radius),
    if radius > 1,

%---- compute norm-factor for parzening 

      nFactor = 1.0; 
		PF(1) = 1.0;
      for line = 2:radius,

        u = (line/(radius+1))*pi;
        v = (sin(u)/u)^2.0;
        PF(line) = v;
        nFactor = nFactor+(2*v*v);
      end;
      nFactor =sqrt(nFactor);

%      -- parzening line 1 

      PF(1) = 1.0/nFactor;

%     -- parzening line 2 to radius      
	for line = 2:radius,
        PF(line) = PF(line)/nFactor;
     end;
     ParInt = sum(PF(2:radius))*2 + PF(1);
PF = PF/ParInt;
else,
   PF = 1;
   end;
   PF= [PF(radius:-1:2) PF]';
   
   
   %----trying to make parzen average a robust est 25.5.3''
   
function[C] = rbpz(S),

sig_m = 1.483*median(abs(S-median(S))); % median absolute deviation 
c_m = 1.5*sig_m;								% Scale
k = find(abs(S) <= c_m);					% assigning Huber weights 
w_m(k) = 1;
q = find(abs(S) > c_m);
w_m(q) = c_m./abs(S(q));					

C = sum(S.*w_m')/sum(w_m);  			% Weighted X spectra

function[TF] = ind_tf(SPMatrix),
j = sqrt(-1);
a = size(SPMatrix);

if length(a) == 4,
    for i = 1:a(1),
  
    XX = squeeze(SPMatrix(i,:,1,1));
    ZX = conj(squeeze(SPMatrix(i,:,2,1)));
    %ZX = (squeeze(SPMatrix(i,:,2,1)));
    YY = squeeze(SPMatrix(i,:,2,2));
    TF(i,:) = ZX./XX; %input/output
    %TF(i,:) = 1./((ZX)./YY); %output/input
    end;
elseif length(a) == 3,

    XX = squeeze(SPMatrix(:,1,1));
    ZX = conj(squeeze(SPMatrix(:,2,1)));
    %ZX = (squeeze(SPMatrix(:,2,1)));
    YY = squeeze(SPMatrix(:,2,2));
    TF = ZX./XX; %input/output
    %TF = 1./((ZX)./YY); %output/input
  
end;

function[Coh] = cohr(SPM),

a = size(SPM);

for i = 1:a(1),
    for j = 1:a(2),
        Coh(i,j) = abs(SPM(i,j,1,2))^2/(SPM(i,j,1,1)*SPM(i,j,2,2));
        if Coh(i,j) > 1,
            Coh(i,j) = 1/Coh(i,j);
        end;
    end;
end;


function[] = disptf(C,Err,freq1,c),
    T = (24*3600)./freq1;
    nfrq = length(C);
   plot(log10(T'),real(C),[c '.-'],'MarkerSize',20);
   hold on;
   plot(log10(T'),imag(C),[c '.-'],'MarkerSize',20);
for i = 1:nfrq,
   tt = real(C(i));
   kk = real(C(i))+Err(i)/2;
   vv = 2*tt-kk;
   h=line(log10([T(i) T(i)]),[kk vv]);
   set(h,'color',c);
   tt = imag(C(i));
   kk = imag(C(i))+Err(i)/2;
   vv = 2*tt-kk;
   h=line(log10([T(i) T(i)]),[kk vv]);
      set(h,'color',c);

end;
